Replication: Young, Old, Conservative, and Bold

This notebook replicates key aspects of Young, Old, Conservative, and Bold by Gârleanu and Panageas published in the Journal of Political Economy in 2015.


The road map for this notebook is:

  1. a brief introduction to the model
  2. a discussion of the equilibrium equations
  3. auxiliary functions and (calibrated) constants
  4. the finite difference methods

Introduction

The model evolves around two types of households, both have Epstein-Zin-Weil preferences. $A$ is bold, takes on more risk and reacts more to a change in the interest rate. On the contrary, $B$ is conservative. Any household has an income that is hump-shaped in its age and death occurs according to a Poisson process with constant intensity, $\pi$, and they have no bequest motives.

Production, $Y_t$, follows a Brownian motion with drift and takes place in a firm. A constant fraction, $\omega$, of production is paid as income to the households, and the rest is distributed as dividends. There is in total one share of this firm and its price is $S_t$. There is also a risk-free bond in zero net supply that has a return $r_t$.

Equilibrium Conditions

The equilibrium is described by the process for the consumption share of type-A households, $X_t$, which is the only state variable:

$$dX_t = \mu_X(X_t) dt + \sigma_X(X_t)dB_t$$

Two coupled ODEs

The above drift and volatility are functions of the solution objects, $g^A$, $g^B$, $\phi^1$, and $\phi^2$. $g^i$ is the consumption-to-wealth ratio and $\phi=\phi^1+\phi^2$ is the ratio of the present value of earnings at birth to aggregate consumption. There are two $\phi^j$ functions because the hump-shape of life-cycle income is parametrized by the sum of two exponential functions.

The equilibrium is described by the following system of coupled, second-order ODEs:

\begin{equation} 0=\frac{\sigma_X^2}{2}\frac{d^2\phi^j}{dX_t^2}+\frac{d\phi^j}{dX_t}\big(\mu_X + \sigma_X(\sigma_Y-\kappa)\big) + \phi^j ( \mu_Y - r - \pi - \delta_j -\sigma_Y \kappa) + B_j \omega \end{equation}

\begin{equation} 0=\frac{\sigma_X^2}{2} M_1^i \Biggm ((M_1^i-1) \biggm(\frac{\frac{dg^i}{dX_t}}{g^i} \biggm)^2+\frac{\frac{d^2g^i}{dX_t^2}}{g^i} \Biggm ) + M_1^i\frac{\frac{dg^i}{dX_t}}{g^i}\big(\mu_X-M_2^i\sigma_X\kappa\big) + \\ \Biggm ( \frac{\kappa^2(X_t)}{2}M_2^i(M_2^i-1)-M_2^i\big(r(X_t)+\pi\big)-M_1^ig^i+\frac{\Xi_3^i}{\gamma^i}\Biggm) \end{equation}

where $\gamma^i$ is relative risk aversion, $\alpha^i = 1 - \frac{1}{\psi^i}$, $\psi^i$ is the elasticity of intertemporal substitution, $\kappa$ is the Sharpe ratio, $\delta^j$ and $B_j$ parametrize the life-cycle income, and $\omega$ is the fraction of production that is paid as income. All other objects are constants that are defined below.

Mind: Equation (A.22) in the paper, which is equation (3) here, has a typo at the second derivative of $g^i$.

The boundary conditions are, for $X_t = 0$ or $X_t = 1$: $$0 = \frac{d\phi^j}{dX_t} \mu_X + \phi^j\big(\mu_Y - r - \pi - \delta_j - \sigma_Y\kappa\big) + B_j \omega$$

$$0 = M_1^i \frac{1}{g^i}\frac{dg^i}{dX_t}\mu_X + \Big(\frac{\kappa(X_t)^2}{2}M_2^i(M_2^i - 1) - M_2^i(r(X_t) + \pi) - M_1^i g^i + \frac{\Xi_3^i}{\gamma^i}\Big)$$

Reformulation

It turns out after a lot of sleuthing, that formulating the system of equations with $g^i$ is much less stable than defining $p^i = \frac{1}{g^i}$ and translating it. The fractions in equation (2) above and the definitions below are substituted with:

$$\frac{\frac{dg^i}{dX_t}}{g^i} = -\frac{\frac{dp^i}{dX_t}}{p^i}$$ and $$\frac{\frac{d^2g^i}{dX_t^2}}{g^i} = 2\Bigg(\frac{\frac{dp^i}{dX_t}}{p^i}\Bigg)^2-\frac{\frac{d^2p^i}{dX_t^2}}{p^i}$$

Auxiliary Functions

The auxiliary functions that appear in the ODEs are defined the following way:

$$\sigma_{X}(X_{t})=\frac{X_{t}\big(\Gamma(X_{t})-\gamma^{A}\big)}{\frac{\Gamma(X_{t})}{\gamma^{B}}X_{t}(1-X_{t})\Big[\frac{1-\gamma^{A}-\alpha^{A}}{\alpha^{A}}\frac{g^{A'}}{g^{A}}-\frac{1-\gamma^{B}-\alpha^{B}}{\alpha^{B}}\frac{g^{B'}}{g^{B}}\Big]+\gamma^{A}}\sigma_{Y}$$

In [1]:
σ_X(X_t, pA, pAprime, pB, pBprime) = σ_Y * X_t * (Γ(X_t) - γA) / 
    (Γ(X_t)/γB * X_t * (1 - X_t) * 
        ((1 - γA - αA) / αA * (-pAprime / pA) - (1 - γB - αB) / αB * (- pBprime / pB)) 
    + γA);

$$\mu_X(X_t) = X_t\Big[\frac{r(X_t)-\rho}{1-\alpha^A} + n^A(X_t)-\pi-\mu_Y\Big] + v^A\pi\beta^A(X_t)-\sigma_Y \sigma_X(X_t)$$

In [2]:
μ_X(X_t, pA, pAprime, pB, pBprime, ϕ1, ϕ2) =  X_t * 
    ((r(X_t, pA, pAprime, pB, pBprime, ϕ1, ϕ2) - ρ) / (1 - αA) + nA(X_t, pA, pAprime, pB, pBprime) - π - μ_Y) +
    vA * π * β(pA, ϕ1, ϕ2) - σ_Y * σ_X(X_t, pA, pAprime, pB, pBprime);

$$\kappa(X_t) = \Gamma(X_t)\sigma_Y + \sum_i\omega^i(X_t)\Big(\frac{1-\gamma^i-\alpha^i}{\alpha^i}\Big) \frac{g^{i'}}{g^i}\sigma_X(X_t)$$

In [3]:
κ(X_t, pA, pAprime, pB, pBprime) =  Γ(X_t) * σ_Y + 
    ωA(X_t) * (1-γA-αA) / αA * (-pAprime / pA) * σ_X(X_t, pA, pAprime, pB, pBprime) +
    ωB(X_t) * (1-γB-αB) / αB * (-pBprime / pB) * σ_X(X_t, pA, pAprime, pB, pBprime);

$$r(X_{t})=\rho+\frac{1}{\Theta(X_{t})}\left\{ \mu_{Y}-\pi\Big(\sum_{i}v^{i}\beta^{i}(X_{t})-1\Big)\right\} -\frac{1}{\Theta(X_{t})}\sum_{i}X_{t}^{i}n^{i}(X_{t})$$

In [4]:
r(X_t, pA, pAprime, pB, pBprime, ϕ1, ϕ2) = ρ + 
    1 / Θ(X_t) * (μ_Y - π * (vA * β(pA, ϕ1, ϕ2) + vB * β(pB, ϕ1, ϕ2) - 1)) -
    1 / Θ(X_t) * (X_t * nA(X_t, pA, pAprime, pB, pBprime) + (1 - X_t) * nB(X_t, pA, pAprime, pB, pBprime));

Mind: The second line in equation (A.16) is not correct, see this issue on GitHub!

Properly derived the definition of nA and nB is:

$$n^{i}(X_{t})=\frac{2-\alpha^{i}}{2\gamma^{i}(1-\alpha^{i})}\kappa^{2}(X_{t})+\frac{\alpha^{i}+\gamma^{i}-1}{2\gamma^{i}\alpha^{i}}\Bigg(\frac{g^{i'}}{g^{i}}\sigma_{X}(X_{t})\Bigg)^{2}+\frac{\alpha^{i}+\gamma^{i}-1}{\alpha^{i}\gamma^{i}}\Big(\frac{g^{i'}}{g^{i}}\sigma_{X}(X_{t})\Big)\kappa(X_{t}) $$

In [5]:
nA(X_t, pA, pAprime, pB, pBprime) =  (2 - αA) / (2 * γA * (1-αA)) * κ(X_t, pA, pAprime, pB, pBprime)^2 + 
    (αA + γA - 1) / (2 * γA * αA) * ((-pAprime / pA) * σ_X(X_t, pA, pAprime, pB, pBprime))^2 +
    (αA + γA - 1) / (γA * αA) * ((-pAprime / pA) * σ_X(X_t, pA, pAprime, pB, pBprime)) * κ(X_t, pA, pAprime, pB, pBprime)
        

nB(X_t, pA, pAprime, pB, pBprime) =  (2 - αB) / (2 * γB * (1-αB)) * κ(X_t, pA, pAprime, pB, pBprime)^2 + 
    (αB + γB - 1) / (2 * γB * αB) * ((-pBprime / pB) * σ_X(X_t, pA, pAprime, pB, pBprime))^2 +
    (αB + γB - 1) / (γB * αB) * ((-pBprime / pB) * σ_X(X_t, pA, pAprime, pB, pBprime)) * κ(X_t, pA, pAprime, pB, pBprime)
Out[5]:
nB (generic function with 1 method)

$$\beta^i(X_t)=g^i(X_t)\underset{\phi^2(X_t) + \phi^2(X_t)}{\underbrace{\phi(X_t)}}$$

In [6]:
β(p, ϕ1, ϕ2) =  (1 / p) * (ϕ1 + ϕ2);

$$X_t^A = X_t\quad\text{and}\quad X_t^B = 1-X_t$$

$$\Gamma(X_t) = \frac{1}{\sum_i \frac{X_t^i}{\gamma^i}}$$

In [7]:
Γ(X_t) =  1 / (X_t / γA + (1 - X_t) / γB);

$$\Theta(X_t) = \sum_i \frac{X_t^i}{1-\alpha^i}$$

In [8]:
Θ(X_t) =  X_t / (1 - αA) + (1 - X_t) / (1 - αB);

$$\omega^i(X_t) = X_t^i \frac{\Gamma(X_t)}{\gamma^i}$$

In [9]:
ωA(X_t) = X_t * Γ(X_t) / γA
ωB(X_t) = (1 - X_t) * Γ(X_t) / γB;

$$\Delta(X_t) = \sum_i \omega^i(X_t)\frac{\gamma^i + 1}{\gamma^i}$$

In [10]:
Δ(X_t) =  ωA(X_t) * (γA + 1) / γA + ωB(X_t) * (γB + 1) / γB;

Constants

The primitive constants of the model are:

In [11]:
const vA = 0.01; const vB = 1 - vA;
const ρ = 0.001
const δ1 = 0.0525; const δ2 = 0.0611
const ω = 1 - 0.08

const γA = 1.5; const ψA = 0.70; const αA = 1 - 1/ψA
const γB = 10.; const ψB = 0.05; const αB = 1 - 1/ψB
const μ_Y = 0.02; const σ_Y = 0.041
const μ = 0.02; const σ = 0.041
const π = 0.02;

Applying the normalization of the integral on the right-hand side in equation (4) in the paper to one yields:

In [12]:
const B1 = 30.72 / (π / (π + δ1) * 30.72 + π / (π + δ2) * -30.29)
const B2 = -30.29 / (π / (π + δ1) * 30.72 + π / (π + δ2) * -30.29);

The auxiliary constants are defined by:

$M_1^i = - 1 - \frac{\Xi_1^i}{\gamma^i}$, $M_2^i = \frac{\gamma^i - 1}{\gamma^i}$, $\Xi_1^i = -\frac{\alpha^i+\gamma^i-1}{\alpha^i}$, $\Xi_2^i = \frac{\alpha^i}{(1-\alpha^i)(1-\gamma^i)}$, $\Xi_3^i = -\frac{\rho + \pi}{\alpha^i}(1-\gamma^i)$, $\Xi_4^i = - \frac{\alpha^i + \gamma^i - 1}{(1-\alpha^i)(1-\gamma^i)}$

In [13]:
const ΞA_1 = - (αA + γA - 1) / αA; const ΞB_1 = - (αB + γB - 1) / αB
const ΞA_2 = αA / ((1-αA) * (1-γA)); const ΞB_2 = αB / ((1-αB) * (1-γB))
const ΞA_3 = - (ρ + π) / αA * (1 - γA); const ΞB_3 = - (ρ + π) / αB * (1 - γB)
const ΞA_4 = - (αA + γA - 1) / ((1-αA) * (1-γA)); const ΞB_4 = - (αB + γB - 1) / ((1-αB) * (1-γB))
const MA_1 = - 1 - ΞA_1 / γA; const MB_1 = - 1 - ΞB_1 / γB
const MA_2 = (γA - 1) / γA; const MB_2 = (γB - 1) / γB;

Numerical Solution - Finite Differences

The following function carries out finite differencing. It takes the evaluation of a solution function on a grid as input and calculates the first derivative via upwinding and the second derivative via central differences. Hence, f, fx, fxx, and wind are vectors with length N.

In [14]:
function diff!(f, fx, fxx, wind, N)
    for i in 1:N
        #fx via upwinding
        if wind[i] >= 0
            fx[i] = (f[min(i + 1, N)] - f[i]) * (N - 1)
        else
            fx[i] = (f[i] - f[max(i - 1, 1)]) * (N - 1)
        end
        
        #fxx not via upwinding
        fxx[i] = (f[min(i + 1, N)] * (N - 1)^2 + f[max(i - 1, 1)] * (N - 1)^2 - 2 * f[i] * (N - 1)^2)
    end
end
Out[14]:
diff! (generic function with 1 method)

Now, the right-hand side of equations (1) for $i\in\{A,B\}$ and equation (2) for $j\in\{1,2\}$ are implemented. The following function takes a stacked candidate solution of the four solution functions evaluated on a grid as input and returns the residual of those four equations.

It carries out the following steps:

  1. unstack the input into the four functions
  2. calculate derivatives without upwinding
  3. calculate the drift $\mu_X$
  4. calculate derivatives via upwinding
  5. return the residuals
In [15]:
function F(Y)
    #unstack the input into the four solution functions
    N = round(Int64, length(Y)/4)
    x = collect(linspace(0, 1, N))
    ϕ1 = Y[1:N]
    ϕ2 = Y[N+1:2*N]
    pA = Y[2*N+1:3*N]
    pB = Y[3*N+1:4*N]
    
    #setup derivatives
    ϕ1x = similar(ϕ1)
    ϕ1xx = similar(ϕ1)
    ϕ2x = similar(ϕ2)
    ϕ2xx = similar(ϕ2)
    pAx = similar(pA)
    pAxx = similar(pA)
    pBx = similar(pB)
    pBxx = similar(pB)

    #finite differences without upwinding
    diff!(ϕ1, ϕ1x, ϕ1xx, zeros(200), N)
    diff!(ϕ2, ϕ2x, ϕ2xx, zeros(200), N)
    diff!(pA, pAx, pAxx, zeros(200), N)
    diff!(pB, pBx, pBxx, zeros(200), N)

    #update the wind 
    wind = μ_X.(x, pA, pAx, pB, pBx, ϕ1, ϕ2)
    diff!(ϕ1, ϕ1x, ϕ1xx, wind, N)
    diff!(ϕ2, ϕ2x, ϕ2xx, wind, N)
    diff!(pA, pAx, pAxx, wind, N)
    diff!(pB, pBx, pBxx, wind, N)

    #the ODEs, stacked into one single vector
    vcat(
        1/2 * ϕ1xx .* σ_X.(x, pA, pAx, pB, pBx).^2 + 
            ϕ1x .* (μ_X.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) + σ_X.(x, pA, pAx, pB, pBx) .* (σ_Y - κ.(x, pA, pAx, pB, pBx))) +
            ϕ1 .* (μ_Y - r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - π - δ1 - σ_Y .* κ.(x, pA, pAx, pB, pBx)) +
            B1 * ω,
        1/2 * ϕ2xx .* σ_X.(x, pA, pAx, pB, pBx).^2 + 
            ϕ2x .* (μ_X.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) + σ_X.(x, pA, pAx, pB, pBx) .* (σ_Y - κ.(x, pA, pAx, pB, pBx))) +
            ϕ2 .* (μ_Y - r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - π - δ2 - σ_Y .* κ.(x, pA, pAx, pB, pBx)) +
            B2 * ω,
        pA .* (1 ./ pA + 
            ψA * (r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - ρ) + κ.(x, pA, pAx, pB, pBx).^2 * (1 + ψA) / (2 * γA) + (1 - ψA * γA) / (γA * (ψA - 1)) * κ.(x, pA, pAx, pB, pBx) .* (pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx)) - 
            (1 - γA * ψA) / (2 * γA * (ψA - 1)) * (pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx)).^2 - π + pAx ./ pA .* μ_X.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) + 0.5 * pAxx ./ pA .* σ_X.(x, pA, pAx, pB, pBx).^2 + 
            (κ.(x, pA, pAx, pB, pBx) / γA + (1 - γA * ψA) / (γA * (ψA - 1)) * (pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx))) .* (pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx)) - r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - κ.(x, pA, pAx, pB, pBx) .* ((pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx)) + (κ.(x, pA, pAx, pB, pBx) / γA + (1 - γA * ψA) / (γA * (ψA - 1)) * (pAx ./ pA .* σ_X.(x, pA, pAx, pB, pBx))))),
        pB .* (1 ./ pB + 
            ψB * (r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - ρ) + κ.(x, pA, pAx, pB, pBx).^2 * (1 + ψB) / (2 * γB) + (1 - ψB * γB) / (γB * (ψB - 1)) * κ.(x, pA, pAx, pB, pBx) .* (pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx)) - 
            (1 - γB * ψB) / (2 * γB * (ψB - 1)) * (pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx)).^2 - π + pBx ./ pB .* μ_X.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) + 0.5 * pBxx ./ pB .* σ_X.(x, pA, pAx, pB, pBx).^2 + 
            (κ.(x, pA, pAx, pB, pBx) / γB + (1 - γB * ψB) / (γB * (ψB - 1)) * (pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx))) .* (pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx)) - r.(x, pA, pAx, pB, pBx, ϕ1, ϕ2) - κ.(x, pA, pAx, pB, pBx) .* ((pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx)) + (κ.(x, pA, pAx, pB, pBx) / γB + (1 - γB * ψB) / (γB * (ψB - 1)) * (pBx ./ pB .* σ_X.(x, pA, pAx, pB, pBx)))))
    )
end
Out[15]:
F (generic function with 1 method)

Expicit Scheme: Newton-Ralphson

One approach is to solve $F(Y) = 0$ by means of guessing an initial value for $Y$ and finding $Y_{t+1}$ from:

$$0=F(Y_t) + J_F(Y_t)\Big(Y_{t+1} - Y_t\Big)$$

This usually needs a good starting value, though. Gârleanu and Panageas solve the ODE with this approach, but use a good starting value.

The following function starts with an initial guess for $Y$ and iterates on the above equation, by calculating the Jacobian via automatic differentiation:

In [16]:
using ForwardDiff

function ExplicitTimeStepping(;initial_Y = ones(800), verbose = true)

    Yold = initial_Y
    Ynew = zeros(Yold)

    distance = 0.1
    iteration = 0

    while (distance > 1e-10) && (iteration < 20)
        Ynew .= -ForwardDiff.jacobian(F, Yold)^-1 * F(Yold) + Yold

        distance = maximum(abs.(Yold .- Ynew))
        Yold = copy(Ynew)
        iteration += 1
        if verbose && (iteration % 1 == 0)
            println("Iteration $iteration with distance $distance")
        end
    end
    
    return Ynew
end
Out[16]:
ExplicitTimeStepping (generic function with 1 method)
In [17]:
ExplicitTimeStepping();
Iteration 1 with distance 585.7008536045665
Iteration 2 with distance 1266.4905411082393
Iteration 3 with distance 641.7484749982939
Iteration 4 with distance 100472.24092114485
Iteration 5 with distance 100182.679406979
Iteration 6 with distance 1139.330470595215
Iteration 7 with distance 6265.674272423209
Iteration 8 with distance 10018.108691268566
Iteration 9 with distance 265557.34131973993
Iteration 10 with distance 264808.42154261214
Iteration 11 with distance 3.671174914817901e6
Iteration 12 with distance 3.669818011934196e6
Iteration 13 with distance 1.774600105757439e6
Iteration 14 with distance 1.1285732151361147e6
Iteration 15 with distance 71193.0326498084
Iteration 16 with distance 52869.319441309235
Iteration 17 with distance 80988.38562556001
Iteration 18 with distance 214153.48411886377
Iteration 19 with distance 992752.2214351855
Iteration 20 with distance 2.0115618457520984e7

This is highly unstable. The starting value is not good enough!

Implicit Scheme: Time-Stepping

Another approach is to solve $F(Y) = 0$ by means of guessing an initial value for $Y$ and finding $Y_{t+1}$ from:

$$0=F(Y_{t+1}) + \frac{1}{\Delta}\Big(Y_{t+1} - Y_t\Big)$$

where $\Delta$ is the step size.

Matthieu Gomez suggests solving this equation via a Newton-Ralphson method in his package EconPDEs, see his notes.

The proposed algorithm has two nested iterations:

  1. The outer iteration over $t$ only concerns the $Y_t$. The outer iteration:

    • begins with a starting value $Y_t = Y_0$
    • runs the inner iteration and receives $Y_{t+1}$ from it
    • if the inner iteration was successful it:
      • calculates the residual of the equation above and aborts if close to zero
      • increases the step size that is to be used in the next iteration
      • updates $Y_t$ with the solution of the inner iteration.
    • if the inner iteration was successful it:
      • discards $Y_{t+1}$
      • decreases the step size
  2. The inner iteration takes $Y_t$ and the step size as given. The above equation is solved exactly for $Y_t+1$ by a Newton-Ralphson method. That is, it iterates over $i$ and begin with $Y_{t+1}^0 = Y_t$ as starting value:

$$ 0 =F(Y_{t+1}^{i})+\frac{1}{\Delta}\Big(Y_{t+1}^{i}-Y_{t}\Big)+ \Big(J_{F}(Y_{t+1}^{i})-\frac{1}{\Delta}\Big)\Big(Y_{t+1}^{i+1}-Y_{t+1}^{i}\Big) $$

This is the first-order Taylor approximation of the above equation expanded around $Y^i_{t}$ and it can be easily solved for $Y^i_{t+1}$. Mind, the Jacobian is a banded matrix, which could be exploited for higher computational performance.

The following function encloses $Y_t$ and the step size. That is, it returns an anonymous function that knows the values of those two objects. The returned function takes a candidate solution of the inner iteration Yprime and an residual vector as input and calculates the residual:

In [18]:
function residual_wrapper(Y, stepsize = 1.)
    return (residual, Yprime) -> residual .= F(Yprime) - 1/stepsize * (Yprime - Y)
end
Out[18]:
residual_wrapper (generic function with 2 methods)

The following function conducts the two nested iteration. The inner iteration is simply a call of the nlsolve function with method = :newton.

In [19]:
using NLsolve

function ImpicitTimeStepping(;initial_Y = ones(800), initial_stepsize = 1., verbose = true)
    
    Yold = initial_Y
    Ynew = zeros(Yold)
    stepsize = initial_stepsize

    distance = 0.1
    distanceold = 0.1
    iteration = 0
    
    while (distance > 1e-10) && (iteration < 200) && (stepsize  >= 1e-12)
        iteration += 1
        
        result = nlsolve(residual_wrapper(Yold, stepsize), Yold, iterations = 25, 
            autodiff=:forward, method = :newton, ftol = 1e-9)
        Ynew .= result.zero
        
        if any(isnan.(Ynew))
            println("Iteration $iteration bad solution with step size $stepsize")
            stepsize = stepsize / 10
            continue
        end
        
        if !result.f_converged
            println("Iteration $iteration no convergence with step size $stepsize")
            stepsize = stepsize / 10
            continue
        end
        
        distance, distanceold = maximum(abs.(F(Ynew))), distance
        Yold = copy(Ynew)
        if (distance <= distanceold)
            stepsize = stepsize * 10
        end
        if verbose && (iteration % 1 == 0)
            println("Iteration $iteration with distance $distance")
        end
    end
    
    return Ynew
end
Out[19]:
ImpicitTimeStepping (generic function with 1 method)

Despite using a naive starting value of 800 ones, the algorithm performs very well:

In [20]:
@time Y = ImpicitTimeStepping();
Iteration 1 with distance 26.110471827973843
Iteration 2 with distance 24.0298426622421
Iteration 3 with distance 14.042368850806449
Iteration 4 with distance 2.1454259551291948
Iteration 5 with distance 0.05369812370985727
Iteration 6 with distance 0.00017256271840082604
Iteration 7 with distance 9.525484756522928e-8
Iteration 8 with distance 4.725109192804666e-12
 26.406420 seconds (33.76 M allocations: 7.000 GiB, 2.37% gc time)

In comparison to the explicit time-stepping method, starting with a candidate that is less than one percent off from the solution maximally (according to the discretization) fails to converge sometimes:

In [ ]:
# scatter plot↔

Evaluation of the Equilibrium

In [21]:
# plot solution functions↔

Plotly javascript loaded.

To load again call

init_notebook(true)

Out[21]:

The following figure replicates Figure 1 from the paper. The volatility of stock prices is different, which yields a different equity premium. This may be a result of the error in equation (A.16) in the paper.

In [22]:
# plot other functions↔
Out[22]:
In [23]:
# plot the law of motion for X↔
Out[23]:

Simulation

Now, the process for $X_t$ is simulated in order to find its stationary distribution:

In [151]:
# plot long path of X↔
Out[151]:

Table 1 in the paper shows the assets pricing relevant objects, which are calculated here again:

In [159]:
# draw table↔
Out[159]:
ParameterBaseline
1Equity Premium0.03704
2Volatility of Returns0.103854
3Average Interest Rate0.0106553
4Volatility of Interest Rate0.000475422

It appears that the model implied parameters are lower than reported in the paper.